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Abstract. In this paper a general theory for interpolation methods on a rectangular grid is 
introduced. By the use of this theory an efficient B-spline based interpolation method for spectral 
codes is presented. The theory links the order of the interpolation method with its spectral properties. 
In this way many properties like order of continuity, order of convergence and magnitude of errors can 
be explained. Furthermore, a fast implementation of the interpolation methods is given. We show 
that the B-spline based interpolation method has several advantages compared to other methods. 
First, the order of continuity of the interpolated field is higher than for other methods. Second, 
only one FFT is needed whereas e.g. Hermite interpolation needs multiple FFTs for computing 
the derivatives. Third, the interpolation error almost matches the one of Hermite interpolation, a 
property not reached by other methods investigated. 

AMS subject classifications. 65T40, 65D05 

Key words. Interpolation, B-spline, three-dimensional, Hermite, Fourier, spline 

1. Introduction. In recent years many studies on the dynamics of inertial par- 
ticles in turbulence have focussed on the Lagrangian properties, see the review by 
Toschi and Bodenschatz pQ . For particles in turbulence, but also in many other appli- 
cations in fluid mechanics, interpolation methods play a crucial role as fluid velocities, 
rate of strain and other flow quantities are generally not available at the location of 
the particles, while these quantities are needed for the integration of the equations of 
motion of the particles. 

When a particle is small, spherical and rigid its dynamics in non-uniform flow 
is governed by the Maxey- Riley (MR) equation [2]. An elaborate overview of the 
different terms in the MR equation and their numerical implementation can be found 
in the paper by Loth [3] and a historical account was given in a review article by 
Michaelides [J]. The evaluation of the hydrodynamic force exerted on the particles 
requires knowledge of the fluid velocity, its time derivative and gradients at the par- 
ticle positions and turns out to be rather elaborate. First, the Basset history force is 
computationally very expensive. However, a significant reduction of cpu-time can be 
obtained by fitting the diffusive kernel of the Basset history force with exponential 
functions, as recently shown by Van Hinsbcrg et al. [5]. Second, the interpolation 
step itself can be very time consuming and memory demanding. Especially for light 
particles, which have a mass density similar to the fluid density (which is, for example, 
sediment transport in estuaries and phytoplankton in oceans and lakes), most terms in 
the Maxey-Riley equation cannot be ignored and therefore also the first derivatives of 
the fluid velocity are needed [6] . For this reason simulations of light particles are com- 
putationally expensive while simulations of heavy particles are less expensive. In order 
to achieve convergence of the statistical properties (probability distribution functions, 
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correlation functions, multi-particle statistics, particle distributions) many particles 
are needed and this calls for fast and accurate interpolation methods. Therefore, our 
aim is to reduce the computation time for the evaluation of the trajectories of light 
particles substantially and make the algorithm competitive with the fast algorithms 
for the computation of trajectories of heavy particles in turbulence. 

The incompressible Navier-Stokes equations arc used to describe the turbulent 
flow field. In turbulence studies the Navier-Stokes equations are often solved by 
means of dealiased pseudo-spectral codes because of the advantage of exponential 
convergence of the computed flow variables. Therefore, we will focus here on interpo- 
lation methods for spectral codes. 

There are many interpolation methods available [7]. We are interested in those 
interpolation methods which are characterized by the following properties. First, the 
method must be accurate, thus we need a high order of convergence. Second, the 
interpolant must have a high order of continuity C p , with p the order of continuity. 
Third, the method must be fast, i.e. computationally inexpensive. A very simple 
interpolation method is linear interpolation. This method is very fast, but unfortu- 
nately this method is relatively inaccurate and it has a low order of convergence. High 
order of convergence can be reached by employing Lagrange interpolation [8] . This in- 
terpolation method has the drawback that it still has a low order of continuity for the 
interpolant. A solution for this problem was recently found by Lalescu et al. [9] who 
proposed a new spline interpolation method. Here, the interpolant has a higher order 
of continuity, but the order of convergence has decreased. A method that has both a 
high order of convergence and a high order of continuity is Hermite interpolation |10| . 
The major disadvantage of this method is that also the derivatives of the function to 
be interpolated are needed, these are calculated by additional Fast Fourier Transforms 
(FFTs) making this method computationally expensive. A remedy to this is B-spline 
interpolarion which has a high order of convergence and errors comparable with 
the ones of Hermite interpolation. Furthermore, this method has a higher order of 
continuity compared with the other methods mentioned above. Normally, the trans- 
formation to the B-spline basis is an expensive step, but by making use of the spectral 
code it can be executed in Fourier space which makes it inexpensive. By executing 
this step in Fourier space the method can be optimized, resulting in smaller errors. 
We believe that the proposed combination of B-splinc interpolation with a spectral 
code makes the method favorable over other traditional interpolation methods. 

Besides exploring the advantages of B-splinc interpolation we focus on neces- 
sary conditions allowing general 3D-intcrpolations to be efficiently executed as suc- 
cessive ID-interpolations. These conditions also carry over desired properties (like 
order of convergence and order of continuity) from the ID-interpolation to the three- 
dimensional equivalent. Further, we provide a fast, generic algorithm to interpolate 
the function and its derivatives using successive ID-interpolations. 

We provide expressions for the interpolation errors in terms of the Fourier compo- 
nents. For this we use Fourier analysis where the interpolation method is represented 
as a convolution function. By doing this, the errors can be calculated as a function of 
the wave number. This gives insight in the behavior of interpolation, especially which 
components are dominant in the interpolation. 

The present study may also be useful for many other applications. Examples in- 
clude the computation of charged particles in a magnetic field |12j US] , but also digital 
filtering and applications in medical imaging [7J [2] . In the latter case interpolations 
are used to improve the resolution of images. Many efforts have been taken to find 
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interpolation methods with optimum qualities [TJ. Still, it is a very active field of 
research. Besides the optimization of interpolation algorithms (accuracy, efficiency), 
the impact of different interpolation methods on physical phenomena like particle 
transport has been investigated in many studies [151 1161 117] . 

In Section [2] we introduce the general framework and explain some one-dimen- 
sional interpolation methods. In Section [3] the framework is generalized to three- 
dimensional interpolation, and a generic algorithm is proposed for the implementa- 
tion of the interpolation in Section 2] A Fourier analysis of the interpolation operator 
is discussed in Section [5j In Section [5] the Fourier analysis is extended to Hermitc 
interpolation and a proof of minimal errors is given. In Section [7j our B-spline based 
interpolation method is introduced, and is compared with three other methods (in- 
cluding Hermite interpolation) in Section [5] Finally, concluding remarks are given in 
Sectional 

2. Interpolation methods. We present a general framework to discuss the var- 
ious interpolation methods. The goal of any interpolation method is to reconstruct 
the original function as closely as possible. As in many applications also some deriva- 
tives of the original function are needed, we focus on reconstructing them as well. 
We start with one-dimensional (ID) interpolation and subsequently, in Section [31 we 
generalize our framework to the three-dimensional (3D) case. 

Let u(x) be a ID function that needs to be reconstructed with x G [0,1]. In 
practice we only have the values of a on a uniform grid, with grid spacing Ax and 
knots at positions Xj, with < j < (Ax) -1 . After interpolation, the function u is 
obtained which is an approximation of u. Now let I be the interpolation operator, so 
u = X[u]. 

When u has periodic boundary conditions, it can be expressed in a Fourier series 
as follows 

u{x) = ]T u k fa{x), fa{x) = e 2 * ikx , (2.1) 
fcez 

with i the imaginary unit and k the wave number. As the grid spacing is finite, a finite 
number of Fourier modes can be represented by the grid. From now on we consider 
u to have a finite number of Fourier modes, so 

fcmax 

u(x) = ^2 ii k (f>k(x), (2.2) 

where fc max , related to Ax, is the maximum wave number. As we add a finite number 
of continuously differentiable Fourier modes fa we have u £ C°°(0,1), a property 
which can be used when constructing the interpolation method. In principle u could 
be reconstructed at any point by the use of Fourier series, however in practice this 
would be far too time consuming and it is therefore not done, instead an interpolation 
is performed, fa is defined as the intcrpolant of fa, i.e., fa = X [fa]. We restrict 
ourselves to linear interpolation operators, i.e., T\ot\U\ + CH2U2] = ail[ui] + a-iL\a^ 
with ai,Q!2 G C This property can be used to write u(x) as 

fcmax 

u{x) = ^2 u k fa{x). (2.3) 

t:=-fc max 

We focus on reconstructing u by pieccwise polynomial functions of degree N — 1 . 
For each interval [xj,Xj+\) with < j < (Ax) -1 we have 
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Here, the vector a depends on the interval under consideration and a T denotes the 
transpose of a. The degree of the highest polynomial function for which the interpo- 
lation is still exact is denoted by n. In this way we get the restriction n < N — 1. 
We consider the reconstruction of u between the two neighboring grid points xj and 
Xj + i. Without loss of generality we can translate and rescale x so that the interval 
[xj ;, Xj+i] becomes the unit interval [0, 1]. 

For Hermite interpolation the values of u and of its derivatives, up to the order 
N/2 — 1 (N must be even), must coincide with those of u at x = and x = 1, i.e., 
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If the derivatives are known then n = N — 1. When the derivatives are not known 
exactly on the grid they have to be approximated by finite difference methods, as 
done by Lalescu et al. [5]. Unfortunately, this method is less accurate than Hcrmitc 
interpolation and n = N — 2. 

The general framework will be illustrated with cubic Hermite interpolation for 
which N = A. So the interpolation uses the function value and the first derivative 
in the two neighboring grid points to construct the interpolation polynomial. We 
have chosen this method because it is very accurate. Moreover, the second derivative, 
which is a piecewise linear function, gives minimal errors with respect to the L 2 -norm. 
This property is further discussed in Section [5] 

First, the discrete values of u and possible derivatives which are given on the grid, 
are indicated with the vector b. In general we have 



b = fM, 



(2.6) 



where the linear operator f depends on the interpolation method and maps a function 
onto a iV-vector. Second, the coefficients of the monomial basis need to be com- 
puted. Because I and f are linear operators, we can write without loss of generality, 



b T M. 



(2.7) 



Here, M is the matrix that defines the interpolation method. For illustration, f and 
M for cubic Hcrmitc interpolation, arc given by 
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Finally, substituting relation (|2.7[) in (|2.4[) gives 



l[u](x) = u(x) = a T x = b T Mx. 



(2.9) 
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In many applications also derivatives are needed. In order to compute the l-th deriva- 
tive of u, the monomial basis functions should be differentiated I times. In general 
this can be done by multiplying a by the differentiation matrix D I times, so 



a T D\ 



D 
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(2.10) 
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where a^) contains the coefficients for the l-th. derivative, obtaining 

^H(ar) = a« T x = b T MDx = b T M«x, 

dx 



(2.11) 



with M (i) = MD 1 . Note that the matrix D is nilpotent, since T> 1 = for I > N, 
implying that at most N — 1 derivatives can be approximated. 

In conclusion, we presented a framework that is able to describe interpolation 
methods, which can be used to implement the interpolation methods in a straightfor- 
ward way. In Section U it is used to generate fast algorithms for the implementation 
of the method. 

3. 3D interpolation. In this section the ID interpolation methods of Section [2] 
are extended to the 3D case. Now the scalar field u depends on the vector x and a 3D 
interpolation needs to be carried out. Like before the interpolated field is given by u 
and I3 is the 3D equivalent of I, so u = I3 [u] . The 3D field u can be represented by 
a 3D Fourier series like 



^M k </> k (x), 



(3.1) 



where </>k is given by 



0k (x) 



^27rik-x 



0k* O)0fc„ (y)(f>kA z ), 



(x,y,z), (3.2) 



and <pk defined by (|2.1[) . Again we restrict ourselves to linear interpolation operators, 
therefore u can be written as 



U(x) = y^M k 0k(x), 
k 



(3.3) 



with </>k the interpolant of <f>\^, i.e., 0k = I3 [0k]- 

The 3D interpolation for a scalar field is carried out applying three times ID 
interpolations, see Fig. 13.11 The interpolation consists of three steps, in which the 
three spatial directions are interpolated one after each other. The order in which the 
spatial directions are interpolated does not matter. Building the 3D interpolation 
out of ID interpolations saves computing time. It can be done for all interpolation 
methods as long as the following two conditions are met. First, the operator I3 must 
be linear. Second, the following condition must be satisfied 



2^3 [0k] =13 [0fc x 0fc„0fcj = 1- [0fc J 1 [0fc v ] I [0fc J =0fc x 0fc B 0fc 2 , (3-4) 



G 




Fig. 3.1. Graphical description of the 3D Lagrange interpolation, using three steps of ID 
interpolations for the case N = 4. First, N 2 ID interpolations are carried out in the x-direction 
(crosses). Second, N interpolations are carried out in y- direction (dots in the right figure) and from 
these N results finally one interpolated value is derived in z-direction (triangle). 



which is the case for almost all interpolation methods. Property Q3.4p can be used 
to prove that properties of the ID operator I carry over to the 3D operator X 3 , for 
example, the order of convergence and the order of continuity. 

Next, relations (|2.9p and (|2.1ip are extended to the 3D case. Like before, we start 
with storing some values of u (given by the spectral code) and possible derivatives in 
the third-order tensor B. In the same fashion as relation (|2.6|) one gets 



B 



(3.5) 



where one element of tensor B is defined like 



Bi 



f y[ f xMii], 



(3.6) 



f x , f y and f z are similar to operator f but now working in a specified direction. For 
Hermite interpolation they arc given by 



u(l,y,z) 



fy[u] 



( u(x,0,z) \ 
u(x, 1, z) 



f.[«] 



/ u(x,y,0) \ 
u(x,y,l) 



(3.7) 



The interpolation is carried out in a similar way as sketched in Fig. 13.11 Similarly to 
()2.9|) . u(x) can be represented as 



X 3 M(x) = u(x) = Bx 1 (Mx)x 2 (My)x 3 (Mz), 



(3.8) 



where M is still the matrix for ID interpolation, y and z arc defined like x which is 
given by relation (|2.4[) . Further, x„ denotes the n-mode vector product [T5], like 



A = Bx n f, 



-l»n + l-"W 



— y ] Bi 1 ---i Ar fi n , 



(3.9) 



where TV denotes the order of tensor B. In this way tensor A is one order less than 
tensor B. Because we employ three of these n-mode vector products the third-order 
tensor B reduces to a scalar. Furthermore, each of these n-mode vector products 
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corresponds to an interpolation in one direction, see also Fig. 13.11 For a general 
derivative one gets 

U : (x) = BXi (M (i) x) x 2 (m ( %) x 3 (m (?c) z) . (3.10) 



dx l dyWz k 

Note that the matrix M does not necessarily have to be the same for the different 
directions x, y and z. One could choose different interpolation methods when for 
example Chebyshev polynomials are used in one direction. In this case the grid is 
nonuniform in this direction and therefore not all interpolation methods can be used. 

Finally, when the scalar field it(x) becomes a vector field u(x), the three compo- 
nents of u can be interpolated separately. This can be written in short by a fourth 
order tensor B where the last dimension contains the three components of u. In this 
way the equations for the new tensor B remain the same as given above. 

Table 4.1 

Algorithm for interpolation, with an estimate of the computational costs 



Computed variables 


Number of flops Number of flops 




for N=4 



x, y and z 3N 12 

Mx, My and Mz 3N 2 48 

M (1) x, M (1) y and M (1) z 3N(N - 1) 36 

Bxi(Mx) 3iV 3 192 



192 



Bxi(Mx)x 2 (My) 3N 2 48 

Bxi(Mx)x 2 (M (1) yj 3N 2 48 

Bxj (M (1 »xj x 2 (My) 3N 2 48 

Bxi(Mx)x 2 (My)x 3 (Mz) 37V 12 



Bxi(Mx)x 2 (My)x 3 (M (1) z) 37V 12 

Bx 1 (Mx)x 2 (M (1) yjx 3 (Mz) 3 AT 12 

Bxi (M (1) x) x 2 (My)x 3 (Mz) 37V 12 

Total: 6iV 3 + 15N 2 + 12 AT 672" 



4. Implementation. Relations (|3.8jl and (|3.10j) provide a good starting point 
for an efficient implementation of the interpolation. We focus on interpolating a 
3D vector field u(x) and on calculating all its first derivatives (which are needed in 
many applications like the computation of the trajectories of inertial particles). The 
matrices M and M^ 1 -* only need to be computed once, which can be done prior to 
interpolation. Second, the vectors x, y and z have to be computed which only needs 
to be done once for each position of interpolation. In Table FTTI we keep track of all the 
computed quantities. Here, the computational costs for evaluating all the components 
is shown where one flop denotes one multiplication with one addition. We show the 
number of flops for the general case and for A^ = 4. The main idea is to reduce the 
order of the tensors as soon as possible in order to generate an efficient method. 

In order to determine how efficient the algorithm is, one can compare the com- 
putational costs against a lower bound. The lower bound we use is related to the size 
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of B which is 37V 3 for a vector field u. In order to be able to use all the information 
in tensor B, 37V 3 flops are needed. For large TV one finds that the algorithm of Table 
14.11 is only a factor 2 less efficient than this lower bound. 

We also compare our algorithm with the one proposed by Lekien and Marsden 
[19] , which uses Hcrmite interpolation with N = 4. Our algorithm has less restrictions 
and shows a slightly better computational performance (for N = 4). The algorithm 
of Lekien and Marsden consists of two steps. First, they calculate the coefficients 
for the polynomial basis. Second, the values at the desired location are calculated. 
They claim that their method is beneficial when the derivatives are needed or the 
interpolation needs to be done multiple times for the same interval, because only the 
second step needs to be executed multiple times. Our method does not have the 
first step, therefore it has no restrictions, nevertheless the computation of the values 
and the first derivatives is slightly faster than for Lekien and Marsden, even when 
considering only the second step. The total costs of their second step is bounded 
by 127V 3 flops (4 times 3iV 3 flops, for the computation of the values and the first 
derivatives). From Table |4~T1 we can conclude that our method needs less flops for the 
same computations. 

5. Fourier analysis. In this section the interpolation operator 1 is expressed 
in terms of a convolution. In this way properties of the interpolation method like the 
order of continuity of the interpolated field and the magnitude of the errors can be 
shown in the Fourier domain. We start with the interpolation of ID functions and 
subsequently, it can be extended to the 3D case. 

Before we start with the derivation, we rescale the variable x by dividing it by Ax, 
so that the new grid spacing equals unity. From now on we work with the rescaled 
grid where x € [0,m] and m = (Ax) , so Xj = j for < j < m. Furthermore we 
introduce the dimensionless wave number k = kAx and <p K is similarly defined as 
4>k, see (|2.1|) . For Hermite interpolation the derivation is somewhat more complex 
because also the derivatives are used and therefore it is postponed to Section [6l We 
focus on interpolation methods where f[u] contains the values of u at the N nearest 
grid points xj of x with local ordering. Thus bj = u(xj) and Xj is given by 



N 

x- -+J 



1,2, •••,7V, (5.1) 



where [-J denotes the nearest lower integer. The interpolation methods can be de- 
scribed by the matrix M, with elements Mjj, see relation (|2.9|) . This relation can 
also be written as 

N 

u(x) = Cj (x — Xj)u (xj) , (5-2) 



with Xj defined by (|5.ip and where Cj is given by 

\ J elsewhere. 

Relation (|5.2j) can be rewritten by using the sifting property of the delta function, 
like 

"0*0 = c i ( x ~ x i) I u (y)Hv - x j)dy- (5.4) 
3=1 J -°° 
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Fig. 5.1. Sketch of linear interpolation as a convolution. The pins represent delta functions 
with the height equal to its prefactor. On the left side is a visualization in real space and on the 
right side in Fourier space. 



This can be further reformulated by subtracting the argument of the delta function 
from the argument of Cj , as 



u (y) ^2 S (y~ Xj)Cj(x - y)dy 

3=1 

u(y)D(y)C(x-y)dy, (5.5) 
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with C[x) and D(x) given by 

TV 

C(x) =Y,Cj{x), D(x) = ^25(x-i). (5.6) 

3=1 i£Z 

In relation (|5.5[) the delta functions can be replaced by the function D, which is a 
train of delta functions because the functions Cj only have a support of length unity, 
see (|5.3p . Finally, the interpolation can be written like 

u = (uD) * C, (5.7) 

with * denoting the convolution product. Here, the convolution function C depends 
on the interpolation matrix M, see Fig. 15.11 

As a consequence of relation (|5.7[) , if the function C is continuous up to the p-th 
derivative then u is also continuous up to the p-th derivative. Even stronger, the order 
of continuity of the function C is equal to the order of continuity of u. Furthermore, 
by the use of relation (|5.7p exact interpolation can be constructed as well 0- 

In the following of this section we will discuss the interpolation error. Before 
proceeding we need to proof the following theorem. 

Theorem. (e K , e\) 2 = for k^A. Here e K is the error in mode k, e K = <f) K — <p K 
and (-)2 is the inner product related to the L 2 -norm || • H2 defined by 

(f,9h= / f(x)g*(x)dx, \\ff 2 = (f,f) 2 = / f(x)f*(x)dx. (5.8) 



The asterisk (*) denotes complex conjugation. 

Proof. We start with replacing u by <j> K in relation (|5.7|) . i.e., 



<j> K = 1 [0 K ] = {<f> K D) * C. (5.9) 
Second, we take the Fourier transform of 4> K , for some fixed kq, i.e., 



(k) = t[ {4>, D) * C] (k) = (F [<t> K0 ] * T[D]) (k)T[C](k) 

J25(k-(i + K Q ))F[C](i + K ), (5.10) 



m 

ie2 



which results in a train of delta functions with the perfector given by J~[C], sec Fig. 
5.11 and J-[-] denotes the Fourier transform given by 



T{g](k) := / g(x)e- 2 ^ x dx. (5.11) 
J — 00 

For linear interpolation these functions are shown in Fig. 15.11 Next, (e Kl e\)2 can 

be written as (e K ,e A ) 2 = (4> K ,4>\j - ^0 k ,^a^) 2 - ^ <j> K , 2 + {4>k,4>\) 2 - Trivially 

(0k,^a) 2 = for k 7^ A. Furthermore, <j) K consists of a discrete set of Fourier com- 
ponents, see relation (|5.10[) . Using this relation, one can show that no common 



1 Exact interpolation can be accomplished by J-"[C](fc) = 1 for —0.5 < k < 0.5 and zero elsewhere. 
In this way only the original Fourier component is filtered out of the spectrum. Note that in this 
case C has infinite support. 
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Fourier components exist for cf> K and (f>\ or cf>\ for k ^ A. Therefore ($> K) <p\^j = 0, 



>\ ) =0 and 



x > = for k ^ A implying (e K , ca)2 = as claimed. 

2 

Corollary. The orthogonality is important to estimate errors. When the error 
in u is computed as \\u — «|||, it can be rewritten like — u||| = ^ K u 2 K ||e K || 2 ) which 
allows easy and straightforward computation of the errors. 

Next, the error in one Fourier component is calculated. In this derivation we 
make use of the fact that 4>k can be written by a sum of Fourier components, see Fig. 
15.11 and relation (|5.10l) . The relative error in one Fourier component can be written 
as 



Un\\t 



'/s|l2 



1 

m 



(T[C] («) 



F[G\ (k + i) (f> K+i 



(5.12) 



From this expression one can see that the error can be computed directly from J-[C] . 



The same can be done for the error in the l-th derivative; e K 



(i) 



. The idea 



is to take the derivatives of the individual Fourier components which results in 



e (0 


2 




2 



= (T[C] («) - 1) 



E 

i/0 



K + i 



[T[C] (« + *)) 5 



(5.13) 



The extension to the 3D case is rather straightforward and is therefore not re- 
ported here. The basic idea is to create 3D functions by multiplying the ID compo- 
nents, this can be done for all functions and the basic equations remain the same. 

6. Hermite interpolation. In this section we extend the theory of Section [5] 
to Hermite interpolation. We also show some special properties that hold for Hermite 
interpolations. Especially, we examine the case N = 4. For this case the second 
derivative becomes a piecewise linear function. Comparison with the actual second 
derivative shows that this piecewise linear function is optimal with respect to the 
L 2 -norm. 

In order to extend the theory of Section [S] to Hermite interpolation with even 
N the same procedure is followed as in Section [S] Analogous to (|5.2p , u(x) can be 
written as 



1 N/2 

l ^ = J2 E C i' 1 ( x - x i) 7^ ( X A ' 

3=0 1=1 



d u 



(6.1) 



where Cjj and Xj are given by 
C jt i(?-3) = 



EtiMt+^j- 1 for0<x<l 
elsewhere 



L*J +3, 



l e i,2,--- , 

3 G 0, 1. 



N 



(6.2) 
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Again following similar steps as in Section [31 u(x) can be rewritten as 

1 N/2 ;_ x 



/OO J I — 1 
-oo 

W /2 ,oo W-l 

= E / d^ (2/)jD(y)G(x ~ y)d2/ ' (6 - 3) 



where D is given by relation (|5.6[) and C; is given by Cj(a:) = Co,j(a;) + Ci.z(x). In 
short, 5 can be written as 

m ( ,i-i \ 



1=1 



similar to relation (|5.7p . Here one can see that for Hcrmite interpolation multiple 
convolution functions Ci are needed which correspond to the derivatives and the 
function itself. Replacing u by K in (|6.4[) gives 

N/2 

4> K = 1[</> K ) = (cb K D) * i^f' 1 Q. (6.5) 
1=1 



In this way we find a similar expression as relation (|5.9I) , where C has to be replaced 
by J2uli (27riK)' _1 C;. In conclusion, relation (|5 . 1 2[) and (|5.13[) can still be used. 
Property. For the error in the first derivative we have the following property: 



e« l) =0, (6.6) 

/ 2 



where the inner product (•, -)2 is defined on the unit interval, i.e. 



(f,g) 2 = [ f(x)g*(x)dx. (6.7) 
Jo 

Furthermore the error in the l-th derivative, e"\ is given by 

e« = |«(x)-|?(*). (6.8) 
da; da; 

Proof. One can rewrite, part of the interpolation conditions for Hcrmite inter- 
polation (|2.5[) in the following way 

/»i , — /• i i 

5(1) - 5(0) = u(l) - u(0) & / ~^dx= —dx. (6.9) 



dx J Q dx 

Here two interpolation conditions give one new condition which is equivalent to rela- 
tion flfEU). 

Corollary. Property (|6.6p shows that the error in the first derivative does not 
have a constant component. Therefore the constant component is exact with respect 
to the L 2 -norm 

Property. For the error in the second derivative in case of N = 4 we have 

e( 2 \l) 2 =0, (e( 2 \a;) 2 = 0. (6.10) 
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Proof. One can rewrite the interpolation conditions (|2.5p for N = 4 in the 
following way 



u(l) -u(0) -u'(0) =u(l) -u(0) -u'(0) & / / —^dxda = / — ^dxda, 

j J dx J Q J dx 



1 ,a d 2~ ,1 .a 



1 j2~ /-I j2 

d u , I a u 



u'(l) - u'{0) = u'(l) - u'(0) & / -r^dx = / -r^dx. (6.11) 

Jo Jo 

The first relation in (|6.10[) follows immediately from the second condition in (|6.11[) . 
The second relation in (|6.10[) is derived in the following way 

Jq Jq = d^^ 

a J {d^ {x) ~d^ {x) ) dx \ a=0 -J (d^ {a) -d^ {a) ) ada = > 

J\^~ d ^) xdx = °- ^ 

Here, the first step is integration by parts and the second step uses the second relation 
of equation (|6.11l) . 

Corollary. Relation (|6.11|) implies that does not have a constant component, 
neither a linear component. For N = 4 the second derivative is a linear function, 
and this means that there is no better approximation in the L -norm of this second 
derivative with a piecewise linear function. This makes Hermite interpolation very 
interesting as a reference case, because we now have proven that the error is minimal 
for this case. 



7. B-spline interpolation. In this section we start with explaining B-splinc 
interpolation. The idea is to create an as smooth as possible intcrpolant. Later it is 
shown how the pseudo-spectral code can be used to efficiently execute this interpo- 
lation method. Furthermore, the interpolation method is optimized to create small 
errors in the L 2 -norm. We start with giving the B-splinc convolution functions af- 
ter which their matrix representation is given and finally the transformation to the 
B-spline basis functions is derived. 

In a spectral code FFTs are applied to transform data from real space to Fourier 
space and backwards. These FFTs are the most expensive step in the simulation and 
therefore we want to keep the number of FFTs needed minimal. This is the reason 
why Hermite interpolation is not a good option, since extra FFTs are needed for the 
computation of the derivatives. An alternative is B-spline interpolation. 

We require high order of continuity of the interpolant. The highest order of conti- 
nuity that can be obtained for the interpolant with piecewise polynomial functions of 
degree N — 1 is C N ~ 2 . In this way the interpolant still matches the original function 
u(x) at the grid points Xj. Moreover, one can immediately see that n = N — 1, where 
n is the highest degree of a polynomial test function for which the interpolation is still 
exact. This high level of continuity can be achieved by using B-splinc functions [TTj . 
The first four uniform B-splinc basis functions Brm are shown in Fig. 17.11 These 
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Fig. 7.1. First four uniform B-splines functions. 



functions can be generated by means of convolutions in the following way 



f 1 for - 0.5 < x < 0.5, 
1 elsewhere, 

B(2) * ^(i), 

B(n) = -B(jv-i) * -B(i)- (7-1) 



These functions have the property that the iV-th function is of degree N — 1 and 
is C N ~ 2 . Furthermore, the B-spline basis functions have local support of length N. 
The B-spline functions can be seen as convolution functions C introduced in Section 
[5] and have a matrix representation. The relation between the functions B/ N \ an d the 
matrix representation is similar to relation (|5.3[) and (|5.6p , and is given by 



-8(2) = 



TV 
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The matrix representation for the first four B-splinc functions is as follows |20| 



M (1) 
M (2) 



(1), 



1 -1 
1 



M 



(3) 



1 

2! 



M (4) = oT 





-2 




(! 


2 


-:•)■ 









/ 1 


-3 


3 -1\ 


4 





-6 3 


1 


3 


3 -3 


V o 





1 J 



(7.3) 



In general we have [20] 
M< 



N 



N-j 



L W,i,j ~ ( N _ ^N- 

with Q l n given by 



iEc- 1 )""*^ - *^-*) 



N-j 



i,i = 1,2,--- ,N, (7.4) 



i\(n — i)\ 



(7.5) 



We still need to express u(x),x 6 Z in terms of B-spline basis functions and 
thus find the transform from real space to the B-spline basis. Because the inverse 
transform from the B-spline basis to real space is somewhat easier, we start with this 
transformation first. From now on we omit the subindex (N). The coefficients of 
the B-spline basis are called ub(x), and u(x) can be derived from it by the discrete 
convolution *u in the following way, u — ub *d Bp. Here, Bp is given by 



B D (x) -- 

and the discrete convolution is given by 



for x < y 



B(x) 

B(x — m) for x > f ' 



x = 0, 1, • • ■ , m — 1, 



(g * D h) (x) = E g{y) h {{x - y) mod 

y=o 



x = 0, 1, • • • , m — 1. 



(7.6) 



(7.7) 



Next, the inverse Bj^ needs to be determined, where B D l is defined by B D l *£> Bo = 
<5d, with Su the discrete delta function, given by 

. , , f 1 for i = „ , oN 

S D (x) = l s = 0,l,-- - ,m-l. (7.8) 

Using the inverse , we can find ub(x) by the discrete convolution ub(x) = u(x)*b. 
B-\x). 

Using a spectral code, the discrete convolution can be evaluated in Fourier space 
and it reduces to a multiplication by constants c(k). These multiplication constants 
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can be computed beforehand and no convolutions need to be evaluated. In this way 
one gets 



To M (k) = F D [u]{k)F D [Bp 1 ] (k) = 
= c{k)F D [u){k), 
where the discrete Fourier transform Tu is given by 



T D \u\(k) 
To [B d ] (k) 



(7.9) 



—2-Kixk/m 



k = 0, 1, • ■ • , m — 1. 



x=0 



(7.10) 



The values of c(k) can be determined in a straightforward manner as suggested above 
using J-d[Bd] but a more optimal choice for c(k) can be made in the following way. 
We minimize the L 2 -norm of the error and for this we use relation (|5.12p and we 
require 



dc(fc) 



<f) K - c(k)(j> K = 0, 



with k = kAx. This implies 



c(k) 



ngj («) 

Y.^AHB] (« + i)) s 



In three dimensions equation ()7.9[) becomes 



7b[«B](k) = c(kx)c(ky)c(k z )J r D [u](k). 



(7.11) 



(7.12) 



(7.13) 



Concluding, we propose an interpolation method for pseudo-spectral codes where 
the interpolation matrix M is given by (|7.4p . Further a multiplication in Fourier 
space is executed like (|7.13[) where the coefficients can be determined from (|7.12[) . 
The coefficients can be computed beforehand and therefore no extra FFTs arc needed, 
making this method very efficient. 

8. Comparison of the interpolation methods. In this section four different 
interpolation methods arc compared. The criteria we are interested in are the follow- 
ing. First, the method must be fast, which is needed because many interpolations will 
usually be carried out. Second, as we are using a spectral code, exponential conver- 
gence is expected and in order to meet this accuracy the interpolation methods must 
have high order of convergence. Furthermore, as the original function is C°°, the in- 
terpolated function must have a high order of continuity as well. Finally, the method 
must have small overall errors. In this way, also the derivatives of the interpolated 
field are still accurate enough. 

The methods that are investigated are the following. First we have Lagrange in- 
terpolation where a polynomial function of degree TV — 1 passes through TV points [8] . 
Second we have investigated the spline interpolation proposed by Lalescu et al. [9]. 
Third, Hermite interpolation is considered and finally our newly proposed B-spline 
interpolation method is used. In Table I8TTI some properties of the interpolation meth- 
ods are reported. All the interpolation methods use piecewise polynomial functions 
of degree TV — 1 to reconstruct the field. 
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method 


n 




order of 


FFT 


comment 








continuity 






Lagrange interpolation 


N - 


1 





1 


for even N 








-1 


1 


for odd TV 


spline interpolation [9] 


N - 


2 


(N-2)/2 


1 


only even iV 


Hermite interpolation 


N - 


1 


(N - 2)/2 


(iV/2) 3 


only even iV 


B-spline based interpolation 


N - 


1 


7V-2 


1 


all N 



Table 8.1 

Overview of the interpolation methods investigated. For all methods the degree of the polynomial 
function is equal to N — 1. 



In order to estimate errors we use relation (f5T2| to find wave number dependent 
errors. For the four interpolation methods these errors are shown in Fig 18.11 In our 
case k max Ax — -| in order to avoid aliasing during the computation of the nonlinear 
term |21| . If this problem were not present, k max could be increased till k max Ax = g. 
From Fig 18.11 also the order of convergence can be determined and it is found equal 
to n + 1 (the lowest degree of a polynomial function for which the interpolation is not 
exact), in agreement with Table I5TT1 




Fig. 8.1. Relative interpolation error for the Fourier mode, see Eg. (5.12p . For all methods 
N = 4 except for linear interpolation which has N = 2. The subfigure shows the Fourier transform 
of two interpolation kernels (spline and B-spline), where the solid line represents exact interpolation. 



In order to avoid extra FFTs the interpolated field can be differentiated as done in 
relation (|3.10[) and the error can be computed by means of (|5 . 13|) . The interpolation 
errors of the first and second derivative are shown in Fig l8.2l Here linear interpolation 
is executed on the derivatives themselves to give a comparison of how accurate the 
interpolation methods arc. One can see for example that the second derivative is 
still better approximated by Hermite interpolation (with N = 4) than with the linear 
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interpolation executed on the second derivative. 
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Fig. 8.2. Relative interpolation error for the first (left) and second derivative (right). Here the 
linear interpolation is taken on the first and second derivative where the other methods are taken 
on the function itself and then differentiated afterwards. Again all methods are taken with N = 4 
except for linear interpolation which has N = 2. 

When comparing the interpolation methods one can see that all interpolation 
methods have a weak point on one of our criteria except for the B-spline based method. 
The Lagrange interpolation for example is only C° continuous for even N and even 
discontinuous for odd N . Furthermore, the overall error is relatively high compared 
with the other methods. The spline interpolation has already a better order of conti- 
nuity but it has a lower order of convergence. Also the overall error is relatively high 
compared with the other methods. Hermite interpolation on the other hand has an 
excellent overall error, especially for the second derivative, see Section [5] The main 
disadvantage of this method is that multiple FFTs are needed which is very time 
consuming. The B-spline based interpolation does not have this problem. The time 
it takes to execute the multiplication in Fourier space can be neglected compared 
with one FFT. Furthermore, this method reaches a much higher order of continuity 
compared with the other methods. When looking at the overall errors in Fig. I8.ll 
and l8.2l one can see that they almost match the one of Hermite interpolation. This is 
especially interesting for the second derivative because we have proven that there can 
not be a better approximation. Note that this second derivative is still continuous for 
the B-spline interpolation whereas for Hermite interpolation it is not. 

9. Conclusions. We have introduced a general framework for interpolation me- 
thods on a rectangular grid. Making use of this framework an algorithm is proposed 
for fast evaluation of the interpolation in three dimensions. This can easily save 
considerable computing time compared with other algorithms. It is shown that the 
computation time needed for this algorithm is close to a theoretical lower bound. 

A spectral theory about these interpolation methods is presented, with which the 
spectral properties of the interpolation methods can be studied. Here basic properties 
of the interpolation method were shown like the order of continuity and the order of 
convergence. Furthermore, errors can be calculated for all Fourier components and 
also for its derivatives. By the use of this theory a novel B-spline based interpolation 
method is introduced for application in conjunction with spectral codes. 

Finally, the interpolation methods for spectral codes are compared. The B-spline 
based interpolation method has several advantages compared with traditional me- 
thods. The order of continuity of the interpolated field is higher than that of Hermite 
interpolation and the other methods investigated. Second, only one FFT needs to be 
done where Hermite interpolation needs multiple FFTs for computing the derivatives. 
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Third, the interpolation error matches almost the one of Hcrmite interpolations which 
is not reached by the other methods investigated. The proposed B-spline interpolation 
is thus the preferred candidate for particle tracking algorithms applied for turbulent 
flow simulations. 
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